Практическое задание 7

Моделирование самосборки липидного бислоя

Плискин Александр
In [1]:
import nglview

import IPython.display
import ipywidgets
from IPython.display import display,display_svg,SVG,Image,HTML

import __main__
from time import sleep
import pymol

import numpy as np

import matplotlib.pyplot as plt
In [2]:
%matplotlib inline
In [3]:
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/b.top
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/pr.mdp
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp
--2019-12-12 19:34:57--  http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp
Resolving kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)... 93.180.63.127
Connecting to kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)|93.180.63.127|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 12772 (12K)
Saving to: `dppc.itp'

100%[======================================>] 12,772      --.-K/s   in 0s      

2019-12-12 19:34:57 (107 MB/s) - `dppc.itp' saved [12772/12772]

--2019-12-12 19:34:57--  http://kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp
Resolving kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)... 93.180.63.127
Connecting to kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)|93.180.63.127|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 34337 (34K)
Saving to: `lipid.itp'

100%[======================================>] 34,337      --.-K/s   in 0.001s  

2019-12-12 19:34:57 (36.8 MB/s) - `lipid.itp' saved [34337/34337]

--2019-12-12 19:34:57--  http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro
Resolving kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)... 93.180.63.127
Connecting to kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)|93.180.63.127|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 2329 (2.3K)
Saving to: `dppc.gro'

100%[======================================>] 2,329       --.-K/s   in 0s      

2019-12-12 19:34:57 (351 MB/s) - `dppc.gro' saved [2329/2329]

--2019-12-12 19:34:57--  http://kodomo.cmm.msu.ru/~golovin/bilayer/b.top
Resolving kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)... 93.180.63.127
Connecting to kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)|93.180.63.127|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 573
Saving to: `b.top'

100%[======================================>] 573         --.-K/s   in 0s      

2019-12-12 19:34:57 (47.4 MB/s) - `b.top' saved [573/573]

--2019-12-12 19:34:58--  http://kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp
Resolving kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)... 93.180.63.127
Connecting to kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)|93.180.63.127|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1088 (1.1K)
Saving to: `em.mdp'

100%[======================================>] 1,088       --.-K/s   in 0s      

2019-12-12 19:34:58 (173 MB/s) - `em.mdp' saved [1088/1088]

--2019-12-12 19:34:58--  http://kodomo.cmm.msu.ru/~golovin/bilayer/pr.mdp
Resolving kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)... 93.180.63.127
Connecting to kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)|93.180.63.127|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1164 (1.1K)
Saving to: `pr.mdp'

100%[======================================>] 1,164       --.-K/s   in 0s      

2019-12-12 19:34:58 (189 MB/s) - `pr.mdp' saved [1164/1164]

--2019-12-12 19:34:58--  http://kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp
Resolving kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)... 93.180.63.127
Connecting to kodomo.cmm.msu.ru (kodomo.cmm.msu.ru)|93.180.63.127|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1277 (1.2K)
Saving to: `md.mdp'

100%[======================================>] 1,277       --.-K/s   in 0s      

2019-12-12 19:34:58 (144 MB/s) - `md.mdp' saved [1277/1277]

На основе одного липида созадим ячейку с 64 липидами.

In [4]:
%%bash
genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
editconf -f dppc.gro -o dppc.pdb
editconf -f b_64.gro -o b_64.pdb
Read 50 atoms
Volume: 1.5477 nm^3, corresponds to roughly 600 electrons
No velocities found
Read 3200 atoms
Volume: 99.0529 nm^3, corresponds to roughly 44500 electrons
No velocities found
                         :-)  G  R  O  M  A  C  S  (-:

                  Green Red Orange Magenta Azure Cyan Skyblue

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                               :-)  genconf  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f       dppc.gro  Input        Structure file: gro g96 pdb tpr etc.
  -o       b_64.gro  Output       Structure file: gro g96 pdb etc.
-trj       traj.xtc  Input, Opt.  Trajectory: xtc trr trj gro g96 pdb cpt

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    0       Set the nicelevel
-nbox        vector 4 4 4   Number of boxes
-dist        vector 0 0 0   Distance between boxes
-seed        int    0       Random generator seed, if 0 generated from the
                            time
-[no]rot     bool   no      Randomly rotate conformations
-[no]shuffle bool   no      Random shuffling of molecules
-[no]sort    bool   no      Sort molecules on X coord
-block       int    1       Divide the box in blocks on this number of cpus
-nmolat      int    3       Number of atoms per molecule, assumed to start
                            from 0. If you set this wrong, it will screw up
                            your system!
-maxrot      vector 180 180 180  Maximum random rotation
-[no]renumber  bool yes     Renumber residues


gcq#226: "Jesus Not Only Saves, He Also Frequently Makes Backups." (Myron Bradshaw)

                         :-)  G  R  O  M  A  C  S  (-:

                  Green Red Orange Magenta Azure Cyan Skyblue

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                               :-)  editconf  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f       dppc.gro  Input        Structure file: gro g96 pdb tpr etc.
  -n      index.ndx  Input, Opt.  Index file
  -o       dppc.pdb  Output, Opt! Structure file: gro g96 pdb etc.
-mead      mead.pqr  Output, Opt. Coordinate file for MEAD
 -bf      bfact.dat  Input, Opt.  Generic data file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    0       Set the nicelevel
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-[no]ndef    bool   no      Choose output from default index groups
-bt          enum   triclinic  Box type for -box and -d: triclinic, cubic,
                            dodecahedron or octahedron
-box         vector 0 0 0   Box vector lengths (a,b,c)
-angles      vector 90 90 90  Angles between the box vectors (bc,ac,ab)
-d           real   0       Distance between the solute and the box
-[no]c       bool   no      Center molecule in box (implied by -box and -d)
-center      vector 0 0 0   Coordinates of geometrical center
-aligncenter vector 0 0 0   Center of rotation for alignment
-align       vector 0 0 0   Align to target vector
-translate   vector 0 0 0   Translation
-rotate      vector 0 0 0   Rotation around the X, Y and Z axes in degrees
-[no]princ   bool   no      Orient molecule(s) along their principal axes
-scale       vector 1 1 1   Scaling factor
-density     real   1000    Density (g/L) of the output box achieved by
                            scaling
-[no]pbc     bool   no      Remove the periodicity (make molecule whole again)
-resnr       int    -1       Renumber residues starting from resnr
-[no]grasp   bool   no      Store the charge of the atom in the B-factor
                            field and the radius of the atom in the occupancy
                            field
-rvdw        real   0.12    Default Van der Waals radius (in nm) if one can
                            not be found in the database or if no parameters
                            are present in the topology file
-sig56       real   0       Use rmin/2 (minimum in the Van der Waals
                            potential) rather than sigma/2 
-[no]vdwread bool   no      Read the Van der Waals radii from the file
                            vdwradii.dat rather than computing the radii
                            based on the force field
-[no]atom    bool   no      Force B-factor attachment per atom
-[no]legend  bool   no      Make B-factor legend
-label       string A       Add chain label for all residues
-[no]conect  bool   no      Add CONECT records to a .pdb file when written.
                            Can only be done when a topology is present


gcq#226: "Jesus Not Only Saves, He Also Frequently Makes Backups." (Myron Bradshaw)

                         :-)  G  R  O  M  A  C  S  (-:

                  Green Red Orange Magenta Azure Cyan Skyblue

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                               :-)  editconf  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f       b_64.gro  Input        Structure file: gro g96 pdb tpr etc.
  -n      index.ndx  Input, Opt.  Index file
  -o       b_64.pdb  Output, Opt! Structure file: gro g96 pdb etc.
-mead      mead.pqr  Output, Opt. Coordinate file for MEAD
 -bf      bfact.dat  Input, Opt.  Generic data file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    0       Set the nicelevel
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-[no]ndef    bool   no      Choose output from default index groups
-bt          enum   triclinic  Box type for -box and -d: triclinic, cubic,
                            dodecahedron or octahedron
-box         vector 0 0 0   Box vector lengths (a,b,c)
-angles      vector 90 90 90  Angles between the box vectors (bc,ac,ab)
-d           real   0       Distance between the solute and the box
-[no]c       bool   no      Center molecule in box (implied by -box and -d)
-center      vector 0 0 0   Coordinates of geometrical center
-aligncenter vector 0 0 0   Center of rotation for alignment
-align       vector 0 0 0   Align to target vector
-translate   vector 0 0 0   Translation
-rotate      vector 0 0 0   Rotation around the X, Y and Z axes in degrees
-[no]princ   bool   no      Orient molecule(s) along their principal axes
-scale       vector 1 1 1   Scaling factor
-density     real   1000    Density (g/L) of the output box achieved by
                            scaling
-[no]pbc     bool   no      Remove the periodicity (make molecule whole again)
-resnr       int    -1       Renumber residues starting from resnr
-[no]grasp   bool   no      Store the charge of the atom in the B-factor
                            field and the radius of the atom in the occupancy
                            field
-rvdw        real   0.12    Default Van der Waals radius (in nm) if one can
                            not be found in the database or if no parameters
                            are present in the topology file
-sig56       real   0       Use rmin/2 (minimum in the Van der Waals
                            potential) rather than sigma/2 
-[no]vdwread bool   no      Read the Van der Waals radii from the file
                            vdwradii.dat rather than computing the radii
                            based on the force field
-[no]atom    bool   no      Force B-factor attachment per atom
-[no]legend  bool   no      Make B-factor legend
-label       string A       Add chain label for all residues
-[no]conect  bool   no      Add CONECT records to a .pdb file when written.
                            Can only be done when a topology is present


gcq#226: "Jesus Not Only Saves, He Also Frequently Makes Backups." (Myron Bradshaw)

In [5]:
__main__.pymol_argv = ['pymol','-qc']
pymol.finish_launching()
pymol.cmd.reinitialize()
pymol.cmd.load('b_64.pdb')
pymol.cmd.do('png b_64.png')
In [6]:
Image(filename='b_64.png')
Out[6]:

В файле b.top установите правильное количество липидов в системе (меняем 34 на 64)

In [7]:
! sed -i 's/34/64/g' b.top

Сделаем небольшой отступ в ячейке от липидов, что бы добавить примерно 2500 молекул воды.

In [8]:
%%bash
editconf -f b_64.gro -o b_ec -d 0.5

Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты молекул.

In [9]:
%%bash
grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v

Запустив в консоли получил:
F{max}^0 = 4.37970e+05
F
{max}^{1} = 6.16887e+02

Добавим в ячейку молекулы воды типа spc.

In [10]:
%%bash
genbox -cp b_em -p b -cs spc216 -o b_s

Проведём "утряску" воды:

In [11]:
%%bash
grompp -f em -c b_s -p b -o b_empr -maxwarn 1
mdrun -deffnm b_empr -v
In [12]:
%%bash
grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
In [13]:
%%bash
editconf -f b_pr.gro -o b_pr.pdb
editconf -f b_s.gro -o b_s.pdb
In [14]:
__main__.pymol_argv = ['pymol','-qc'] 
pymol.finish_launching()
pymol.cmd.load('b_s.pdb')
pymol.cmd.do('png b_s.png')
In [15]:
Image(filename='b_s.png')
Out[15]:
In [16]:
__main__.pymol_argv = ['pymol','-qc'] 
pymol.finish_launching()
pymol.cmd.load('b_pr.pdb')
pymol.cmd.do('png b_pr.png')
In [17]:
Image(filename='b_pr.png')
Out[17]:

Заметно, что вода более равномерно и плотно расположилась внутри ячейки. Сама ячейка стала "беспорядочна"

In [18]:
%%bash
grompp -f md -c b_pr -p b -o b_md -maxwarn 1

Анализ результатов вычислений на кластере

In [19]:
%%bash
echo 2 | trjconv -f b_md.gpu.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol
In [20]:
import io
from base64 import b64encode
video = io.open('mo.mp4', 'r+b').read()
video_encoded = b64encode(video)
html_code = '<video controls alt="PyMol Movie" src="data:video/mp4;base64,{}" type="video/mp4">'.format(video_encoded)
display(HTML(html_code))

На 59 изображении отчетливо видно разделение на два слоя.
MODEL59 соответствует t=29000.00000

In [21]:
Image(filename='img59.png')
Out[21]:
In [22]:
%%bash
echo 2 | g_traj -f b_md.gpu.xtc -s b_md.tpr -ob box_1.xvg
In [23]:
import pandas as pd

coord = pd.read_csv('box_1.xvg', header = None, skiprows = 24, sep='\t')
coord = coord[[1,2,3,4]]
coord.columns = ['t', 'x', 'y', 'z']
coord.head(10)
Out[23]:
t x y z
0 0 6.26000 4.44300 5.77800
1 25 6.22935 4.42124 5.81437
2 50 6.20596 4.40464 5.85363
3 75 6.17243 4.38085 5.91878
4 100 6.13236 4.35240 5.96941
5 125 6.13492 4.35422 5.98395
6 150 6.13508 4.35434 6.02683
7 175 6.12250 4.34541 6.02864
8 200 6.12901 4.35003 6.02500
9 225 6.11738 4.34177 6.03904

Площадь считаем как произведение по осям Y и Z и нормируем на 32 (число липидов в слое)

In [24]:
plt.figure(figsize=(15,8))
lip_sizes = np.array(coord.y)*np.array(coord.z)/32.
plt.plot(np.array(coord.t), lip_sizes);

Судя по графику со временем площадь липида уменьшается. В процессе моделирования структура компактно собирается, приходя в оптимальное состояние.

Посмотрим на площадь доступных гидрофильной и гидрофобной поверхностей.

In [25]:
%%bash
echo 2 | g_sas -f b_md.gpu.xtc -s b_md.tpr -o sas_b.xvg -dt 100
In [26]:
hydro = pd.read_csv('sas_b.xvg', header = None, skiprows = 22, sep='\s+')
hydro = hydro[[0,1,2]]
hydro.columns = ['t', 'Hydrophobic', 'Hydrophilic']
hydro.head(10)
Out[26]:
t Hydrophobic Hydrophilic
0 0 206.248 135.565
1 100 184.583 141.929
2 200 179.992 154.704
3 300 157.336 168.297
4 400 148.188 165.633
5 500 138.280 164.318
6 600 132.765 168.017
7 700 137.785 166.972
8 800 130.916 170.924
9 900 127.910 169.047
In [27]:
plt.figure(figsize=(20,8))
plt.plot(np.array(hydro.t), hydro.Hydrophobic, label='Hydrophobic', color='red')
plt.plot(np.array(hydro.t), hydro.Hydrophilic, label='Hydrophilic', color='blue')
plt.legend()
Out[27]:
<matplotlib.legend.Legend at 0x7fe7808208d0>

Площадь доступной гидрофобной поверхности стремительно снижается в первые ~5000 шагов, что соответствует увиденному выше образованию гидрофобного ядра. Это логично: упаковка гидрофобной поверхности в ядро существенно снижает энергию системы. В дальнейшем обе площади поверхностей продолжают снижаться, хотя и менее интенсивно: ядро обращается в бислой, в котором ещё меньше места для контактов с растворителем.

Считаем меру порядка бислоя в начале и в конце моделирования.

Скачиваем индекс файл

In [28]:
! wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx
In [29]:
%%bash
g_order -s b_md -f b_md.gpu.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
g_order -s b_md -f b_md.gpu.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X
                         :-)  G  R  O  M  A  C  S  (-:

                          GROtesk MACabre and Sinister

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                               :-)  g_order  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f   b_md.gpu.xtc  Input        Trajectory: xtc trr trj gro g96 pdb cpt
  -n        sn1.ndx  Input        Index file
 -nr      index.ndx  Input        Index file
  -s       b_md.tpr  Input        Run input file: tpr tpb tpa
  -o  ord_start.xvg  Output       xvgr/xmgr file
 -od     deuter.xvg  Output       xvgr/xmgr file
 -ob      eiwit.pdb  Output       Protein data bank file
 -os     sliced.xvg  Output       xvgr/xmgr file
 -Sg     sg-ang.xvg  Output, Opt. xvgr/xmgr file
 -Sk    sk-dist.xvg  Output, Opt. xvgr/xmgr file
-Sgsl sg-ang-slice.xvg  Output, Opt. xvgr/xmgr file
-Sksl sk-dist-slice.xvg  Output, Opt. xvgr/xmgr file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    19      Set the nicelevel
-b           time   0       First frame (ps) to read from trajectory
-e           time   5000    Last frame (ps) to read from trajectory
-dt          time   0       Only use frame when t MOD dt = first time (ps)
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-xvg         enum   xmgrace  xvg plot formatting: xmgrace, xmgr or none
-d           enum   x       Direction of the normal on the membrane: z, x or y
-sl          int    1       Calculate order parameter as function of box
                            length, dividing the box in #nr slices.
-[no]szonly  bool   no      Only give Sz element of order tensor. (axis can
                            be specified with -d)
-[no]unsat   bool   no      Calculate order parameters for unsaturated
                            carbons. Note that this cannot be mixed with
                            normal order parameters.
-[no]permolecule bool   no      Compute per-molecule Scd order parameters
-[no]radial  bool   no      Compute a radial membrane normal
-[no]calcdist  bool no      Compute distance from a reference (currently
                            defined only for radial and permolecule)

Taking x axis as normal to the membrane
Reading file b_md.tpr, VERSION 4.5.5 (single precision)
Using following groups: 
Groupname: C34 First atomname: C34 First atomnr 33
Groupname: C36 First atomname: C36 First atomnr 35
Groupname: C37 First atomname: C37 First atomnr 36
Groupname: C38 First atomname: C38 First atomnr 37
Groupname: C39 First atomname: C39 First atomnr 38
Groupname: C40 First atomname: C40 First atomnr 39
Groupname: C41 First atomname: C41 First atomnr 40
Groupname: C42 First atomname: C42 First atomnr 41
Groupname: C43 First atomname: C43 First atomnr 42
Groupname: C44 First atomname: C44 First atomnr 43
Groupname: C45 First atomname: C45 First atomnr 44
Groupname: C46 First atomname: C46 First atomnr 45
Groupname: C47 First atomname: C47 First atomnr 46
Groupname: C48 First atomname: C48 First atomnr 47
Groupname: C49 First atomname: C49 First atomnr 48
Groupname: C50 First atomname: C50 First atomnr 49

Reading frame       0 time    0.000   Number of elements in first group: 64
Last frame        200 time 5000.000   

Read trajectory. Printing parameters to file
Atom 1 Tensor: x=-0.00721784 , y=-0.0117279, z=0.0189457
Atom 2 Tensor: x=-0.0170186 , y=-0.00653596, z=0.0235546
Atom 3 Tensor: x=-0.00360557 , y=-0.0130777, z=0.0166833
Atom 4 Tensor: x=-0.00482021 , y=-0.00504665, z=0.00986689
Atom 5 Tensor: x=-0.00630494 , y=-0.00518794, z=0.0114929
Atom 6 Tensor: x=-0.00496223 , y=-0.00843269, z=0.013395
Atom 7 Tensor: x=-0.00207883 , y=-0.00591445, z=0.00799332
Atom 8 Tensor: x=-0.00653175 , y=0.00305707, z=0.00347471
Atom 9 Tensor: x=-0.00764073 , y=0.00344101, z=0.00419976
Atom 10 Tensor: x=-0.0117189 , y=0.0023166, z=0.00940237
Atom 11 Tensor: x=-0.00572851 , y=0.00779041, z=-0.00206186
Atom 12 Tensor: x=-0.00421947 , y=0.00626815, z=-0.00204864
Atom 13 Tensor: x=-0.00230188 , y=0.00412894, z=-0.00182703
Atom 14 Tensor: x=-0.00370037 , y=0.00205719, z=0.00164321

Back Off! I just backed up ord_start.xvg to ./#ord_start.xvg.5#

Back Off! I just backed up deuter.xvg to ./#deuter.xvg.11#

gcq#225: "My Ass May Be Dumb, But I Ain't No Dumbass." (Jackie Brown)

                         :-)  G  R  O  M  A  C  S  (-:

                  Green Red Orange Magenta Azure Cyan Skyblue

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                               :-)  g_order  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f   b_md.gpu.xtc  Input        Trajectory: xtc trr trj gro g96 pdb cpt
  -n        sn1.ndx  Input        Index file
 -nr      index.ndx  Input        Index file
  -s       b_md.tpr  Input        Run input file: tpr tpb tpa
  -o    ord_end.xvg  Output       xvgr/xmgr file
 -od     deuter.xvg  Output       xvgr/xmgr file
 -ob      eiwit.pdb  Output       Protein data bank file
 -os     sliced.xvg  Output       xvgr/xmgr file
 -Sg     sg-ang.xvg  Output, Opt. xvgr/xmgr file
 -Sk    sk-dist.xvg  Output, Opt. xvgr/xmgr file
-Sgsl sg-ang-slice.xvg  Output, Opt. xvgr/xmgr file
-Sksl sk-dist-slice.xvg  Output, Opt. xvgr/xmgr file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    19      Set the nicelevel
-b           time   45000   First frame (ps) to read from trajectory
-e           time   0       Last frame (ps) to read from trajectory
-dt          time   0       Only use frame when t MOD dt = first time (ps)
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-xvg         enum   xmgrace  xvg plot formatting: xmgrace, xmgr or none
-d           enum   x       Direction of the normal on the membrane: z, x or y
-sl          int    1       Calculate order parameter as function of box
                            length, dividing the box in #nr slices.
-[no]szonly  bool   no      Only give Sz element of order tensor. (axis can
                            be specified with -d)
-[no]unsat   bool   no      Calculate order parameters for unsaturated
                            carbons. Note that this cannot be mixed with
                            normal order parameters.
-[no]permolecule bool   no      Compute per-molecule Scd order parameters
-[no]radial  bool   no      Compute a radial membrane normal
-[no]calcdist  bool no      Compute distance from a reference (currently
                            defined only for radial and permolecule)

Taking x axis as normal to the membrane
Reading file b_md.tpr, VERSION 4.5.5 (single precision)
Using following groups: 
Groupname: C34 First atomname: C34 First atomnr 33
Groupname: C36 First atomname: C36 First atomnr 35
Groupname: C37 First atomname: C37 First atomnr 36
Groupname: C38 First atomname: C38 First atomnr 37
Groupname: C39 First atomname: C39 First atomnr 38
Groupname: C40 First atomname: C40 First atomnr 39
Groupname: C41 First atomname: C41 First atomnr 40
Groupname: C42 First atomname: C42 First atomnr 41
Groupname: C43 First atomname: C43 First atomnr 42
Groupname: C44 First atomname: C44 First atomnr 43
Groupname: C45 First atomname: C45 First atomnr 44
Groupname: C46 First atomname: C46 First atomnr 45
Groupname: C47 First atomname: C47 First atomnr 46
Groupname: C48 First atomname: C48 First atomnr 47
Groupname: C49 First atomname: C49 First atomnr 48
Groupname: C50 First atomname: C50 First atomnr 49

Reading frame       0 time 45000.000   Number of elements in first group: 64
Reading frame     190 time 49750.000   

Read trajectory. Printing parameters to file
Atom 1 Tensor: x=0.0435801 , y=-0.00455514, z=-0.0390249
Atom 2 Tensor: x=0.0284894 , y=0.00199356, z=-0.0304829
Atom 3 Tensor: x=0.0220051 , y=0.005912, z=-0.027917
Atom 4 Tensor: x=-0.00375097 , y=0.00509763, z=-0.00134663
Atom 5 Tensor: x=0.00558598 , y=-0.00480467, z=-0.00078128
Atom 6 Tensor: x=-0.000900556 , y=-0.0127421, z=0.0136427
Atom 7 Tensor: x=0.00337893 , y=-0.0146854, z=0.0113065
Atom 8 Tensor: x=-0.0120317 , y=-0.0093706, z=0.0214024
Atom 9 Tensor: x=-0.00964038 , y=-0.0101485, z=0.0197889
Atom 10 Tensor: x=-0.0105375 , y=-0.00752532, z=0.0180629
Atom 11 Tensor: x=-0.0159844 , y=-0.000336015, z=0.0163205
Atom 12 Tensor: x=-0.0114039 , y=-0.00232148, z=0.0137254
Atom 13 Tensor: x=-0.0123366 , y=-0.001215, z=0.0135517
Atom 14 Tensor: x=-0.0079866 , y=-0.00766345, z=0.0156501

Back Off! I just backed up ord_end.xvg to ./#ord_end.xvg.6#

Back Off! I just backed up deuter.xvg to ./#deuter.xvg.12#

gcq#89: "This Doesn't Suck, It's a Black Hole !" (K.A. Feenstra)

In [30]:
start = pd.read_csv('ord_start.xvg', header = None, skiprows = 12, sep='\s+')
end = pd.read_csv('ord_end.xvg', header = None, skiprows = 12, sep='\s+')
In [31]:
plt.figure(figsize=(17,8))
plt.plot(np.array(start[0]), start[3], label='Start', color='red');
plt.plot(np.array(end[0]), end[3], label='End', color='green');
plt.legend();
Out[31]:

К концу моделирования билипидный слой существенно упорядоченнее, чем вскоре после начала.